/*==============================================================================
FILE NAME: Figure_A6.do
CREATED: 13 June 2025
==============================================================================*/

// Panel A
use "$processed_data/Air_Panel_week.dta", clear

drop if never_air_inv == 1

egen t = group(year week)
xtset RN_id t

// Generate variables for post-incident outcomes
forv h = 0/12 {
	gen p_air_inv_`h' = f`h'.p_air_inv - l1.p_air_inv
}

// Generate pre-treatment ΔP(investigation) variables for 2–12 weeks before incident
forv h = 2/12 {
	gen p_air_inv_neg`h' = l`h'.p_air_inv - l1.p_air_inv
}

// Create fixed effect identifiers
egen RN_year = group(RN_id year)
gen quarter = 1 if inlist(month, 1, 2, 3)
replace quarter = 2 if inlist(month, 4, 5, 6)
replace quarter = 3 if inlist(month, 7, 8, 9)
replace quarter = 4 if month > 9

egen RN_quarter = group(RN_id year quarter)
egen RN_month   = group(RN_id year month)

// Keep necessary variables only
keep p_air_inv_* p_air_incident RN_year t RN_id 

// Initialize containers for coefficients and confidence intervals
cap drop b u d se Years Zero
gen Years = _n - 13 if _n <= 25
gen Zero = 0 if _n <= 25
gen b = 0
gen se = 0
gen u = 0
gen d = 0                             

// Run event study regressions: post-treatment weeks
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
    reghdfe p_air_inv_`h' p_air_incident, ///
        absorb(RN_year t) cluster(RN_id)
    unique t if e(sample)
    replace b = _b[p_air_incident] if Years == `h'
    replace se = _se[p_air_incident] if Years == `h'
    replace u = (_b[p_air_incident] + 1.96 * _se[p_air_incident]) if Years == `h'
    replace d = (_b[p_air_incident] - 1.96 * _se[p_air_incident]) if Years == `h'
}

// Run event study regressions: placebo (pre-treatment) weeks
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
    reghdfe p_air_inv_neg`h' p_air_incident, ///
        absorb(RN_year t) cluster(RN_id)
    unique t if e(sample)
    replace b = _b[p_air_incident] if Years == -`h'
    replace se = _se[p_air_incident] if Years == -`h'
    replace u = (_b[p_air_incident] + 1.96 * _se[p_air_incident]) if Years == -`h'
    replace d = (_b[p_air_incident] - 1.96 * _se[p_air_incident]) if Years == -`h'
}

// Prepare and plot the graph
preserve
keep if Years != .
keep b u d se Years Zero
graph set window fontface "Times New Roman"
twoway ///
    (rarea u d Years, col(navy) fint(inten20) lwidth(0) lpattern(solid)) ///
    (line b Years, lcolor(navy) lpattern(solid) lwidth(medium)) ///
    (line Zero Years, lcolor(black)), ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(0(0.1)0.6, labsize(vlarge)) ///
    legend(off) ///
    ytitle("{&Delta} P(Investigation)", size(vlarge)) ///
    xtitle("Week", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)

graph export "$figures/Figure_A6_Panel_A.pdf", replace
restore

// Panel B
use "$processed_data/Air_Panel_week.dta", clear

drop if never_air_inv == 1
egen t = group(year week)
xtset RN_id t

// Re-generate ΔP variables
forv h = 0/12 {
	gen p_air_inv_`h' = f`h'.p_air_inv - l1.p_air_inv
}
forv h = 2/12 {
	gen p_air_inv_neg`h' = l`h'.p_air_inv - l1.p_air_inv
}

// Fixed effects: RN-by-quarter
egen RN_year = group(RN_id year)

gen quarter = 1 if inlist(month, 1, 2, 3)
replace quarter = 2 if inlist(month, 4, 5, 6)
replace quarter = 3 if inlist(month, 7, 8, 9)
replace quarter = 4 if month > 9

egen RN_quarter = group(RN_id year quarter)
egen RN_month   = group(RN_id year month)

// Keep necessary variables
keep p_air_inv_* p_air_incident RN_quarter t RN_id 

// Initialize storage
cap drop b u d se Years Zero
gen Years = _n - 13 if _n <= 25
gen Zero = 0 if _n <= 25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

// Regressions for post-treatment
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
    reghdfe p_air_inv_`h' p_air_incident, ///
        absorb(RN_quarter t) cluster(RN_id)
    unique t if e(sample)
    replace b = _b[p_air_incident] if Years == `h'
    replace se = _se[p_air_incident] if Years == `h'
    replace u = (_b[p_air_incident] + 1.96 * _se[p_air_incident]) if Years == `h'
    replace d = (_b[p_air_incident] - 1.96 * _se[p_air_incident]) if Years == `h'
}

// Regressions for pre-treatment 
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
    reghdfe p_air_inv_neg`h' p_air_incident, ///
        absorb(RN_quarter t) cluster(RN_id)
    unique t if e(sample)
    replace b = _b[p_air_incident] if Years == -`h'
    replace se = _se[p_air_incident] if Years == -`h'
    replace u = (_b[p_air_incident] + 1.96 * _se[p_air_incident]) if Years == -`h'
    replace d = (_b[p_air_incident] - 1.96 * _se[p_air_incident]) if Years == -`h'
}

// Plot Panel B
keep if Years != .
keep b u d se Years Zero
graph set window fontface "Times New Roman"
twoway ///
    (rarea u d Years, col(navy) fint(inten20) lwidth(0) lpattern(solid)) ///
    (line b Years, lcolor(navy) lpattern(solid) lwidth(medium)) ///
    (line Zero Years, lcolor(black)), ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(0(0.1)0.6, labsize(vlarge)) ///
    legend(off) ///
    ytitle("{&Delta} P(Investigation)", size(vlarge)) ///
    xtitle("Week", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)

// Save figure
graph export "$figures/Figure_A6_Panel_B.pdf", replace
